/**
 * @file   channels.hpp
 * @author Kimmo Luoma <kimmo@PHOTON3>
 * @date   Mon Sep 10 15:12:48 2012
 * 
 * @brief  Class for representing jump channels. Keys form a sequence.
 * 
 * 
 */
#ifndef CHANNELS_HPP
#define CHANNELS_HPP
#include "rate.hpp"
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>

#include <map>
#include <utility>
using namespace boost::numeric::ublas;

class channels{
private:
  int dim_of_H;
  int number_of_channels;
  std::map<int,matrix<double> > ops;
  std::map<int,rate> D_rates;
  std::map<int,rate> L_rates;
public:
  /** 
   * Default constructor
   * 
   */
  channels();
  
  /** 
   * Copy constructor
   * 
   * @param chan 
   */
  channels(const channels& chan);
  
  /** 
   * Copy assignment
   * 
   * @param chan 
   * 
   * @return object type channels
   */
  channels& operator=(const channels& chan);
  
  /** 
   * Destructor
   * 
   */
  ~channels();
  
  /** 
   * Set operator \f$ C_k \f$ to the channels. If operator with key k does not yet
   * exist in channels then it is added. If operator with key k
   * exists then it the matrix representation is replaced with new one
   * 
   * @param C_k, Matrix representation of the operator.
   * @param k, key of the operator. 
   */
  void set_op(matrix<double>  C_k, int k);
  
  /** 
   * Set new operator C_k to the channels with largest key,
   * k > k' for all k' already existing in channels.
   *  
   * @param C_k,   matrix representation of the operator.
   */
  void set_op(matrix<double> C_k);
  
  /** 
   * Set decay rate \f$ \Delta_k(t) \f$ for channels. 
   * If rate with key k does not yet
   * exist in channels then it is added. If rate with key ii
   * exists then the rate is replaced with new one
   * 
   * @param D_k rate
   * @param k key for rate
   */
  void set_decay_rate(rate D_k,int k);
  
  /** 
   * Set new decay rate \f$ \Delta_k(t) \f$ to the channels with largest key,
   * k>k', for all k' existing already in channels.
   *    
   * @param D_k
   */
  void set_decay_rate(rate D_k);

  /** 
   * Set Lamb shift rate \f$ S_k(t) \f$ for channels. 
   * If rate with key k does not yet
   * exist in channels then it is added. If rate with key ii
   * exists then the rate is replaced with new one
   * 
   * @param S_k rate
   * @param k key for rate
   */
  void set_lamb_rate(rate S_k,int k);
  
  /** 
   * Set new Lamb shift rate \f$ S_k(t) \f$ to the channels with largest key,
   * k>k', for all k' existing already in channels.
   *    
   * @param S_k
   */
  void set_lamb_rate(rate S_k);
  
  /** 
   * Get the number of channels in the system
   * 
   * 
   * @return number of channels <int>
   */
  int get_number_of_channels();
  
  /** 
   * Get decay rate with key k.
   * 
   * @param k key <int>
   * 
   * @return object of type rate
   */
  rate get_decay_rate(int k);

  /** 
   * Get Lamb shift rate with key k.
   * 
   * @param k key <int>
   * 
   * @return object of type rate
   */
  rate get_lamb_rate(int k);
  
  /** 
   * Get operator with key ii.
   * 
   * @param k Key of the rate
   * 
   * @return rate class.
   */
  matrix<double> get_op(int k);
  
  /** 
   * Return jump operator (time dependent).
   * 
   * 
   * @return \f$ \Delta_k C_k \f$ 
   */
  matrix<double> operator()(int k,double t);
  /** 
   * Return the effect of Lamb Shifts into the Hamiltonian at time t.
   * 
   * 
   * @return \f$ \sum_k S_k(t)C_k^\dagger C_k \f$ 
   */
  matrix<double> operator()(double t);
  /** 
   * Return lamb shift term of coherent evolution
   * 
   * @param t 
   * 
   * @return \f$ \sum_k S_k(t)C_k^\dagger C_k \f$ 
   */
  matrix<double> get_lamb_shift(double t);
  
};
#endif
